Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 February 2008 (MN I^TeX style file vl.4) 



Thermal conduction in cosmological SPH simulations 



> 



Martin Jubelgas^^, Volker Springel^ and Klaus Dolag^ 

^ Max- Planck- Institut fiir Astrophysik, Karl-Schwarzschild-Strafie 1, 85740 Garching bet Miinchen, Germany 
- - - , ^ Dipartimento di Astronomia, Universita di Padova, Viccolo dell' Osservatorio 5, 1-35122 Padova, Italy 

o 

^ , 2 February 2008 

^ . ABSTRACT 

Thermal conduction in the intracluster medium has been proposed as a possible heat- 
I ing mechanism for offsetting central cooling losses in rich clusters of galaxies. However, 

because of the coupled non-linear dynamics of gas subject to radiative cooling and 
thermal conduction, cosmological hydrodynamical simulations are required to reliably 
predict the effects of heat conduction on structure formation. In this study, we intro- 
' duce a new formalism to model conduction in a diffuse ionised plasma using smoothed 

I particle hydrodynamics (SPH), and we implement it in the parallel TreePM/SPH-code 

, GADGET-2. We consider only isotropic conduction and assume that magnetic suppres- 

sion can be described in terms of an effective conductivity, taken as a fixed fraction 
I of the temperature-dependent Spitzer rate. We also account for saturation effects in 

. low-density gas. Our formulation manifestly conserves thermal energy even for indi- 

' vidual and adaptive timesteps, and is stable in the presence of small-scale temperature 

fH ^ noise. This allows us to evolve the thermal diffusion equation with an explicit time 

integration scheme along with the ordinary hydrodynamics. We use a series of simple 
test problems to demonstrate the robustness and accuracy of our method. We then 
apply our code to spherically symmetric realizations of clusters, constructed under the 
assumptions of hydrostatic equilibrium and a local balance between conduction and 
^ I radiative cooling. While we confirm that conduction can efficiently suppress cooling 

flows for an extended period of time in these isolated systems, we do not find a similarly 
strong effect in a first set of clusters formed in self-consistent cosmological simulations. 
However, their temperature profiles are significantly altered by conduction, as is the 
X-ray luminosity. 
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1 INTRODUCTION 

Clusters of galaxies provide a unique laboratory to study 
structure formation and the material content of the Uni- 
verse, because they are not only the largest virialized sys- 
tems but are also believed to contain a fair mixture of cos- 
mic matter. Among the many interesting aspects of cluster 
physics, their X-ray emission takes a particularly prominent 
role. It provides direct information on the thermodynamic 
state of the diffuse intracluster gas, which makes up for most 
of the baryons in clusters. 

While the bulk properties of this gas are well under- 
stood in terms of the canonical ACDM model for structure 
formation, there are a number of discrepancies between ob- 
servations and the results of present hydrodynamical sim- 
ulations. For example, a long standing problem is to un- 
derstand in detail the scaling relations of observed clusters, 
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which deviate significantly from simple self-similar predic- 
tions. In particular, poor clusters of galaxies seem to contain 
gas of higher entropy in their centres than expected (Pon- 
man et al., 1999; Lloyd-Davies et al., 2000). This has been 
interpreted either to be evidence for an entropy injection due 
to non-gravitational processes (Loewenstein, 2000; Wu et al., 
1999; Metzler & Evrard, 1994), or as a sign of the selective 
removal of low-entropy gas by gas cooling (Voit et al., 2002; 
Wu & Xue, 2002). Both processes combined could influence 
the thermodynamic properties of the ICM in a complex in- 
terplay (Tornatore et al., 2003; Borgani et al., 2003). 

Another interesting problem occurs for the radial tem- 
perature profiles of clusters. Most observed clusters show a 
nearly isothermal temperature profile, often with a smooth 
decline in their central parts (Allen et al., 2001; Johnstone 
et al., 2002; Ettori et al., 2002). Nearly isothermal profiles 
are also obtained in adiabatic simulations of cluster forma- 
tion (e.g. Frenk et al., 1999). However, clusters in simula- 
tions that include dissipation typically show temperature 
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profiles that increase towards the centre (e.g. Lewis et al., 
2000), quite different from what is observed. 

Perhaps the biggest puzzle is that spectroscopic X-ray 
observations of the centres of clusters of galaxies have re- 
vealed little evidence for cooling of substantial amounts of 
gas out of the intracluster medium (e.g. David et al., 2001), 
even though this would be expected based on their bolo- 
metric X-ray luminosity alone (Fabian, 1994). The apparent 
absence of strong cooling flows in clusters hence indicates the 
presence of some heating source for the central intracluster 
medium. Among the proposed sources axe AGN, buoyant ra- 
dio bubbles (Churazov et al., 2001; Enfflin & Heinz, 2001), 
feedback processes from star formation (Bower et al., 2001; 
Menci & Cavaliere, 2000), or acoustic waves (Fujita et al., 
2003). 

Recently, Narayan & Medvedev (2001) have proposed 
that thermal conduction may play an important role for the 
cooling processes in clusters. The highly ionised hot plasma 
making up the ICM in rich clusters of galaxies should be effi- 
cient in transporting thermal energy, unless heat diffusion is 
inhibited by magnetic fields. If conduction is efficient, then 
cooling losses in the central part could be offset by a con- 
ductive heat flow from hotter outer parts of clusters, which 
forms the basis of the conduction idea. 

Indeed, using simple hydrostatic cluster models where 
cooling and conductive heating are assumed to be locally 
in equilibrium, Zakamska & Narayan (2003, ZN henceforth) 
have shown that the central temperature profiles of a num- 
ber of clusters can be well reproduced in models with con- 
duction (see also Fabian et al., 2002; Voigt et al., 2002; 
Briiggen, 2003; Voigt & Fabian, 2003). The required con- 
ductivities are typically sub-Spitzer (Medvedev et al., 2003). 
This suggests that thermal conduction may play an impor- 
tant role for the thermodynamic properties of the ICM. 

On the other hand, it has been frequently argued (Chan- 
(iran & Cowley, 1998; Malyshkin & Kulsrud, 2001) that 
magnetic fields in clusters most likely suppress the effective 
conductivity to values well below the Spitzer value for an 
unmagnetised gas. Note that rotation measurements show 
that magnetic fields do exist in clusters (e.g. Vogt & Enfilin, 
2003). However, little is known about the small-scale field 
configuration, so that there is room for models with chaoti- 
cally tangled magnetic fields (Narayan & Medvedev, 2001), 
which may leave a substantial fraction of the Spitzer con- 
ductivity intact. The survival of sharp temperature gradi- 
ents along cold fronts, as observed by Chandra in several 
clusters (Markevitch et al., 2000; Vikhhnin et al., 2001; Et- 
tori & Fabian, 2000), may require an ordered magnetic field 
to suppress conduction. For a more complete review of the 
effects of magnetic fields on galactic clusters, sec Carilli & 
Taylor (2002) and references therein. 

It is clearly of substantial interest to understand in de- 
tail the effects conduction may have on the formation and 
structure of galaxy clusters. In particular, it is far from clear 
whether the temperature profile required for a local balance 
between cooling and conduction can naturally arise during 
hierarchical formation of clusters in the ACDM cosmology. 
This question is best addressed with cosmological hydrody- 
namical simulations that fully account for the coupled non- 
hnear dynamics of the gas subject to radiative cooling and 
thermal conduction. 

In this paper, we hence develop a new numerical imple- 



mentation of conduction and include it in a modern TreeSPH 
code for structure formation. Smoothed particle hydrody- 
namics (Lucy, 1977; Monaghan & Lattanzio, 1985; Mon- 
aghan, 1992) is a powerful numerical tool to investigate gas 
dynamics, which has found widespread application in as- 
trophysics. In cosmological simulations, gas densities vary 
over many orders of magnitude, and can change rapidly 
as function of position and time. Classical mesh-based Eu- 
clidean approaches to hydrodynamics have difficulty adjust- 
ing to this high dynamic range unless sophisticated adaptive 
mesh refinements (AMR) methods are used. In contrast, the 
Lagrangian approach of SPH, thanks to its matter- tracing 
nature, guarantees good spatial resolution in high-density 
regions, while spending little computational time on low- 
density regions of space, where a coarser spatial resolution 
is usually sufficient. 

Another useful aspect of SPH is that the equations that 
describe physical processes can often be translated into an 
SPH form rather intuitively. Below we discuss this in detail 
for the conduction equation, highlighting in particular how 
a number of practical problems with respect to stability can 
be overcome. The final formulation we propose is robust and 
explicitly conserves thermal energy, even when individual 
timesteps for each particle are used. 

After examining idealised test problems to validate our 
implementation of conduction, we apply our code to real- 
izations of the static cluster models of ZN, investigating in 
particular, to what degree conduction may balance cooling 
in these clusters, and for how long this approximate equilib- 
rium can be maintained. We also discuss results for a first 
set of cosmological simulations of cluster formation. We here 
compare simulations that follow only adiabatic gas dynam- 
ics, or only cooling and star formation, with corresponding 
ones that also include conduction. As we will see, thermal 
conduction can lead to a substantial modification of the final 
thermodynamic properties of rich clusters. 

The outline of this paper is as follows. In Section 2, we 
discuss the basics of the conduction equation, and our nu- 
merical approach for discretising it in SPH. In Section 3, 
we show first test runs of conducting slabs, which we also 
use to illustrate various issues of numerical stability. In Sec- 
tion 4, we consider isolated clusters, constructed with an 
initial equilibrium model, while in Section 5 we present a 
comparison of results obtained in cosmological simulations 
of cluster formation. Finally, we summarise and discuss our 
findings in Section 6. 



2 THERMAL CONDUCTION IN SPH 
2.1 The conduction equation 

Heat conduction is a transport process for thermal energy, 
driven by temperature gradients in the conducting medium. 
Provided the mean free path of particles is small compared 
to the scale length of the temperature variation, the local 
heat flux can be described by 

j = -kVT, (1) 

where T{r) gives the temperature field, and k is the heat 
conduction coefficient, which may depend on local proper- 
ties of the medium. For example, in the case of an astro- 
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physical plasma, we encounter a strong dependence of k on 
the temperature itself. 

The rate of temperature change induced by conduction 
can simply be obtained from conservation of energy, viz. 



dw 

'd7 



-V-j, 



(2) 



where u is the thermal energy per unit mass, and p de- 
notes gas density. Eliminating the heat flux, this can also be 
written directly in terms of u, giving the heat conduction 
equation in the form 



du 
dt 



-V-(kVT). 



(3) 



Note that the temperature is typically simply proportional 
to u, unless the mean particle weight changes in the rele- 
vant temperature regime, for example as a result of a phase 
transition from neutral gas to ionised plasma. In the thin 
astrophysical plasmas we are interested in, the strong tem- 
perature dependence of the Spitzer conductivity (see below) 
suppresses conduction in low-temperature gas heavily. For 
all practical purposes we can set u = kT/^-y — l)/it] = c„T, 
where /u = 0.588 mp is the mean molecular weight of a fully 
ionised gas with the primordial mix of helium and hydrogen, 
and Cv is the heat capacity per unit mass. 



2.2 Spitzer Conductivity 

Spitzer (1962) derived the classical result for the heat con- 
ductivity due to electrons in an ionised plasma. It is given 
by 



1.31 UeXek 



1/2 



(4) 



where Tie is the electron density, and Ae the electron mean 
free path. Interestingly, the product UeXe depends only on 
the electron temperature Te, 



AeTle = 



33/'(fcTe)^ 



(5) 



provided we neglect the very weak logarithmic dependence 
of the Coulomb logarithm A on electron density and tem- 
perature, which is a good approximation for clusters. We 
will set InA = 37.8, appropriate for the plasma in clusters 
of galaxies (Sarazin, 1988). The Spitzer conductivity then 
shows only a strong temperature dependence, Ksp oc T®''^, 
and has the value 



keT erg 
lOkeVy cmskeV 



(6) 



Note that the presence of magnetic fields can in princi- 
ple strongly alter the conductivity. Depending on the field 
configuration, it can be suppressed in certain directions, or 
even in all directions in cases of certain tangled fields. The 
field configuration in clusters is not well understood, and it is 
currently debated to what degree magnetic fields suppress 
conduction. We will assume that the modification of the 
conductivity can be expressed in terms of an effective con- 
ductivity, which we parameterise as a fraction of the Spitzer 
conductivity. 

Even in the absence of magnetic fields, the Spitzer con- 
ductivity can not be expected to apply down to arbitrarily 



low plasma densities. Eventually, the scale length of the tem- 
perature gradient will become comparable or smaller than 
the electron mean free path, at which point the heat flux will 
saturate, with no further increase when the temperature gra/- 
dicnt is increased (Cowie & McKee, 1977). This maximum 
heat flux jsat is given by 

1/2 



2knT 

TTTTle 



(7) 



In order to have a smooth transition between the Spitzer 
regime and the saturated regime, we limit the conductive 
heat fiux by defining an effective conductivity (Sarazin, 
1988) in the form 

= l+4?Ae//x - 

Here It = T/\'VT\ is the characteristic length-scale of the 
temperature gradient. 

2.3 SPH formulation of conduction 

At first sight, equation (3) appears to be comparatively easy 
to solve numerically. After all, the time evolution generated 
by the diffusion equation smoothes out initial temperature 
variations, suggesting that it should be quite 'forgiving' to 
noise in the discretisation scheme, which should simply also 
be smoothed out. 

In practice, however, there are two problems that make 
it surprisingly difficult to obtain stable and robust im- 
plementations of the conduction equation in cosmological 
codes. The first has to do with the presence of second 
derivatives in equation (3), which in standard SPH kernel- 
interpolants can be noisy and sensitive to particle disorder. 
The second is that an explicit time integration method can 
easily lead to an unstable integration if large local gradients 
arise due to noise. We will discuss our approaches to solve 
these two problems in turn. 

A simple discretisation of the conduction equation in 
SPH can be obtained by first estimating the heat flux for 
each particle applying standard kernel interpolation meth- 
ods, and then estimating the divergence in a second step. 
However, this method has been shown to be quite sensitive 
to particle disorder (Brookshaw, 1985), which can be traced 
to the effective double-differentiation of the SPH-kernel. In 
addition, this method has the practical disadvantage that an 
intermediate result, the heat flux vectors, need to be com- 
puted and stored in a separate SPH-loop. 

It is hence advantageous to use a simpler SPH discreti- 
sation of the Laplace operator, which should ideally involve 
only first order derivatives of the smoothing kernel. Such a 
discretisation has been proposed before (Brookshaw, 1985; 
Monaghan, 1992), and we here give a brief derivation of it 
in three dimensions. 

For a well-behaved field ^ (a;), we can consider a Taylor- 
series approximation for Y{xj) in the proximity of Y{xi), 
e.g. 



Y{xj) - Y{xi) 



VY 



{Xj -Xi) + 



1 d^Y 



2 dxgdxk 



Oixj-XiY 



(9) 
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Neglecting terms of third and higher orders, we multiply this 
with 



{xj — Xi)'ViW{xj — Xi) 



Xi\ 



(10) 



where W{x) — W{\x\) is the SPH smoothing kernel. Note 
that we choose this kernel to be spherically symmetric and 
normalised to unity. The expression in equation (10) is well 
behaved for Xj — > Xi under these conditions. Introducing 

the notations Xij — xj — Xi and Wij = W{xj — Xi), we 
integrate over all Xj and note that 



I- 



XiiViWi. 



0, 



(a^ij)s (3 



ij )k- 



(11) 



(12) 



So the term linear in VF drops out, and the terms 
involving off-diagonal elements of the Hesse matrix of Y 
vanish, so that the sum over the second order term simply 
reduces to V^Y. We hence end up with 



-2 



/ 



Y{xj 



^^Xij-ViWijd^Xj. 



(13) 



This analytical approximation of the Laplacian can now 
be easily translated into an SPH kernel interpolant. To this 
end, we can replace the integral by a sum over all particles 
indexed by j, and substitute the volume element d^Xj by its 
discrete SPH analogue nij/pj. The values of the field Y{x) 
at the particle coordinates can be either taken as the value of 
the intrinsic particle property that is evolved, Y{xi) = Vj, or 
as a kernel interpolant of these values, Y(xi) = (Yi), where 
for example 

{Y,) =J2Yj^W{xij). (14) 
j 



We then end up with a discrete SPH approximation of the 
Laplace operator in the form 



„ \ ^ m, Yi 



Yi 



XijViWi- 



(15) 



We now consider how this can be applied to the thermal 
conduction problem, where the conductivity may also show 
a spatial variation. Using the identity 

V(KVr) = i [V^(/tT) -TV^ + kV^T] , (16) 

wo can use our result from equation (15) to write down a 
discretised form of equation (3): 



dui _ 
'dt ~ ^ 



rrij {tij + n^) {Tj - T,) _ 
PiPi kijP 



(17) 



This form is antisymmetric in the particles i and j', and 
the energy exchange is always balanced on a pairwise basis, 
i.e. conservation of thermal energy is manifest. Also, it is 

easy to sec that the total entropy can only increase, and 
that heat always flows from higher to lower temperature. 

The conductivities Kj and Hj in equation (17) are effec- 
tively arithmetically averaged. Cleary & Monaghan (1999) 
proposed to make the replacement 

I'^i ^" ^^j 2Kihij 



l^i ~\~ l^j 



(18) 



They showed that this ensures a continuous heat flux even 
in cases when the heat conductivity exhibits a disconti- 
nuity, as for example along the interface between differ- 
ent phases. It is clear then that this modification should 
also behave better when the conductivity changes extremely 
rapidly on small scales, as it can happen for example in 
ICM gas when cool particles get into direct contact with 
comparatively hot neighbours. Indeed, we found this sym- 
metrisation to give numerically more robust behaviour, par- 
ticularly in simulations that in addition to heat conduc- 
tion also followed radiative cooling. Note that since we have 
min(Ki,Kj) < 2KiKj / {Ki + Kj) < 2mm{Ki,Kj), the Cleary & 
Monaghan average stays always close to the smaller of the 
two conductivities involved, to within a factor of two. 



2.4 Numerical implementation details 

We have implemented heat conduction in a new version 
of the massively paraUel TreeSPH-code GADGET (Springel 
et al., 2001b), which is a general purpose code for cosmolog- 
ical structure formation. Unlike earlier public releases of the 
code, the present version, GADGET-2, uses the 'entropy for- 
mulation' of SPH proposed by Springel & Hernquist (2002) 
which conserves both energy and entropy (in regions with- 
out shocks) for fully adaptive smoothing lengths. In this 
formulation, an entropic function 



A=(7-l)- 



.7-1 



(19) 



is evolved as independent variable for each particle, instead 

of the thermal energy per unit mass. Note that A = A{s) 
is only a function of the thermodynamic entropy s per unit 
mass. 

For consistency with this formalism, we need to express 
the heat conduction equation in terms of entropy. This is 
easily accomplished in an isochoric approximation, noting 
that the temperature can be expressed asT = ix/ks Ap'^~^ , 
where /i is the mean molecular weight. This results in 



dAj 
dt 



kB pi-' ^ PiPi \pY~^ 



Ai \ XijViWi- 
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Pi 



(20) 

One problematic aspect of the heat conduction equa- 
tion is that small-scale numerical noise in the temperature 
field can generate comparatively large heat flows, simply be- 
cause this noise can involve small-scale gradients of sizable 
magnitude. Since we are using an explicit time integration 
scheme for the hydrodynamical evolution, this immediately 
raises the danger of instabilities in the integration. Unless 
extremely small timesteps or an implicit integration scheme 
are used, the energy exchange between two particles due 
to a small-scale temperature difference can become so large 
that the explicit time integration "overshoots" , thereby po- 
tentially reversing the sign of the temperature difference be- 
tween the two particles in conductive contact. This is not 
only incorrect, but makes it possible for the temperature 
differences to grow quickly in an oscillatory fashion, causing 
an instable behaviour of the integration. 

We have found that a good method to avoid this prob- 
lem is to use a kernel interpolant for the temperature field 
(or entropy field) in the discretisation of the heat conduc- 
tion equation (20), instead of the individual particle tem- 
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perature values themselves. The interpolant represents a 
smoothed version of the noisy sampling of the temperature 
field provided by the particle values, so that on the scale 
of the SPH smoothing length, small-scale noise in the heat 
flux is strongly suppressed. Particles will still try to equili- 
brate their temperatures even within the smoothing radius 
of each particle, but this will happen at a damped rate. Heat- 
conduction due to temperature gradients on larger scales is 
unaffected however. In a later section of this work, we will 
explicitly demonstrate how this improves the stability of the 
time integration, particularly when individual and adaptive 
timesteps are used. 

For definitcncss, the interpolant we use is a smoothed 
version Ai of the entropy, defined by 



(21) 



We then replace Ai and/or Aj on the right-hand-side of 
equation (20) with the intorpolants Ai and Aj. Note that 
the weighting by p'"'^"'' ensures that we obtain a value of 
A that corresponds to a smoothed temperature field, as re- 
quired since conduction is driven by gradients in temper- 
ature and not entropy. However, since the density values 
need to be already known to evaluate this interpolant, this 
unfortunately requires an additional SPH loop, which causes 
quite a bit of computational overhead. This can in princi- 
ple be avoided if the weighting with densities in equation 
(21) is dropped, which according to our tests appears to be 
sufficiently accurate in most situations. 

Note that although we use the SPH interpolant for the 
entropy values in equation (20), we still compute the values 
for the particle conductivity k based on the intrinsic particle 
temperature and not on its smoothed counterpart. 

While the above formulation manifestly conserves ther- 
mal energy, this property may get lost when individual and 
adaptive time steps are used, where for a given system step 
only a subset of the particles is evolved. There can then be 
particle pairs where only one particle is active in the current 
system step, while the other is not, so that a 'one-sided' 
heat conduction may occur that causes a (brief) violation 
of energy conservation during the step. While most of the 
resulting energy imbalance will only be a temporary fluc- 
tuation that will be compensated as soon as the 'inactive' 
particle in the pair is evolved, the very strong temperature 
dependence of the conductivity may produce sizable errors 
in this situation, particularly when coarse timestepping is 
used. We therefore decided to implement an explicitly con- 
servative scheme for the heat exchange, even when adaptive 
and individual timesteps are used. 

To this end, we define a pairwise exchange of heat en- 
ergy as 



2/i mifnjK.ij 
ks piPj 



A, \ x.jV^Wi^ 
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P 



(22) 



A simple translation into a finite difference scheme for the 
time evolution would then be 



iu'i = ruiUi + Ati ^ Eij, 



(23) 



where u'i and m can also be expressed in terms of the corre- 
sponding entropy values A'^ and Ai. However, for individual 
and variable timesteps, only a subfraction of particles will 



be 'active' in the current system timestep. These particles 
have individual timesteps Ati, while the 'inactive' particles 
can be formally assigned Ati = for the step. Equation (23) 
then clearly does not guarantee detailed energy conservation 
in each system step. 

We recover this property in the following way. We up- 
date the energy of particles according to 



miu'i = rriiUi + i ^ Atj{5i 

jk 



■ SikjEjk 



(24) 



in each system step. In practice, the double-sum on the right 
hand side can be simply computed using the usual SPH loop 
for active particles. Each interacting pair found in the neigh- 
bour search for active particle i is simply used to change the 
thermal energy of particle i by AtiEij/2, and that of the 
neighbouring particle j by —AtiEij/2 as well. Note that if 
particle j is active, it carries out a neighbour search itself 
and will find i itself, so that the total energy change of i due 
to the presence of j is given by {Ati + Atj)Eij/2. Energy 
is conserved by construction in this scheme, independent of 
the values if individual timesteps of particles. If all parti- 
cles have equal steps, equation (24) is identical to the form 
of (23). As before, in our SPH implementation we convert 
the final thermal energy change to a corresponding entropy 
change. These entropy changes are applied instantaneously 
for all particles at the end of one system step, so that even 
particles that are inactive at the current system step may 
have their entropy changed by conduction. The resulting 
low order of the time integration scheme is sufficient for the 
diffusion equation. 

If conduction is so strong that the timescale of conduc- 
tive redistribution of thermal energy is short compared to 
relevant dynamical timescales of the gas, the usual hydrody- 
riamic timestep selected by the code based on the Courant- 
critorion may become too coarse to follow conduction accu- 
rately. We therefore introduce an additional timestep cri- 
terion in simulations that follow thermal conduction. To 
this end, we limit the maximum allowed timestep to some 
prescribed fraction of the conduction time scale, j4/|yicond|, 
namely 



Atcond = ax 



\Ao 



(25) 



where a is a dimensionless accuracy parameter. In our sim- 
ulations, we typically employed a value of a = 0.25, which 
has provided good enough accuracy at a moderate cost of 
CPU time. 

Finally, in order to account for a possible limitation 
of conduction by saturation, we compute the gradient of 
the (smoothed) temperature field in the usual SPH fashion, 
and then use it to compute It and the saturation-limited 
conductivity based on equation (8). 



3 ILLUSTRATIVE TEST PROBLEMS 

In this section, we verify our numerical implementation of 
thermal conduction with a rmmber of simple test problems 
that have known analytic solutions. We will also investigate 
the robustness of our formulation with respect to particle 
disorder, and (initial) noise in the temperature field. 
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Figure 1. Time evolution of tlie temperature profile of two slabs of solid material, brought in contaet with each other at i = along 
the z = 50 cm plane, with an initial difference of thermal energy per unit mass of Am = 1000 erg g~^. Crosses mark numerical results, 
and the solid line is the analytic solution of the heat conduction equation. 



3.1 Conduction in a one-dimensional problem 

We first consider two solid slabs with different initial temper- 
atures, brought into contact with each other at time t = Q. 
The slabs were realized as a lattice of SPH particles, with di- 
mensions 100 X 10 X 10, and an equidistant particle-spacing 
of 1 cm. To avoid perturbations of the ID-symmetry, the 
simulation volume was taken to be periodic along the two 
short axes, while kept non-periodic along the long axis, al- 
lowing us to study surface effects, if present, on the sides of 
the slabs that are not in contact. Note that the test was car- 
ried out with our 3D code. All particle velocities were set to 
zero initially, and kept at this value to mimic a solid body. 
The thermal conductivity weis set to a value corresponding 
to a = t^/(cvp) = 1 cm^s^^ throughout both materials, 
independent of temperature. 

In Figure 1, we compare the time evolution of our nu- 
merical results for this test with the analytical solution of 
the same problem, which is given by 

«(«,t) = Ko + ^erf (^1^), (26) 

where Zm gives the position of the initial difference of size 
A« in thermal energy, and wo is the mean thermal energy. 
We see that the numerical solution tracks the analytic result 
very nicely. 

We have also repeated the test for a particle configura- 
tion corresponding to a 'glass' (White, 1996), with equally 
good results. For a Poisson distribution of particles, we noted 
however a small reduction in the speed of conduction when a 
small number of neighbours of 32 is used. Apparently, here 
the large density fluctuations due to the Poisson process 
combined with the smallness of the number of neighbours 
leads to somewhat poor coupling between the particles. This 



effect goes away for larger numbers of neighbours, as ex- 
pected. Note however that in practice, due to the pressure 
forces in a gas, the typical configuration of tracer particles 
is much more akin to a glass than to a Poisson distribution. 

Next, we examine how robust our formulation is with 
respect to small-scale noise in the temperature field. To this 
end, we repeat the above test using the glass configuration, 
but wc perturb the initial thermal energies randomly with 
fluctuations at an rms- level of cr = 80 erg g"'^. Further, we 
increase the maximum timestep allowed for particles to = 
1.0 s. 

We discussed previously that we can either use the par- 
ticle values of the temperatures (or entropy) in the right- 
hand-side of equation (22), or kernel interpolants thereof. 
Here we compare the following different choices with each 
other: 

(A) Basic formulation: Use particle values Ai and Aj . 

(B) Smoothed formulation: Use Ai and Aj. 

(C) Mixed formulation: Use Ai and Aj . 

The mixed formulation (C) may at first seem problematic, 
because its pair- wise antisymmetry is not manifest. How- 
ever, since we use equation (24) to exchange heat between 
particles, conservation of thermal energy is ensured also in 
this case. When all particles have equal timesteps, formula- 
tion (C) would correspond to an arithmetic average of (A) 
and (B). 

In Figure 2, we compare the result obtained for these 
three formulations after a simulation time of 5 sec. The ar- 
tificially perturbed initial conditions are shown in the top 
left panel. Interestingly, while all three different formula- 
tions are able to recover the analytic solution in the mean, 
they show qualitatively different behaviour with respect to 
the imprinted noise. When the intrinsic particle values for 
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Figure 2. Comparison of tlie numerical stability of the SPH discretisation scheme when different formulations for the conduction equation 

(22) arc used. The top left panel shows the initial conditions for the simple conduction problem considered in Fig. 1, but perturbed 
with artificial Gaussian noise in the temperature field. The top right panel shows the evolved state after 5 sec when formulation (A) is 
employed, where the individual particle temperatures are used directly. The lower left panel compares this with formulation (B), where 
the temperatures on the right hand-side of equation (22) are taken to be kernel interpolants of the temperature field. Finally, the bottom 
right panel gives the result for formulation (C), which represents a mixed scheme that uses both the individual temperatures of particles, 
and the smoothed temperature field. This scheme proves to be the most robust against local noise, which is quickly damped away. 



the temperatures are used (top right panel) , very large pair- 
wise gradients on small scales occur that induce large heat 
exchanges. As a result, the particles oscillate around the 
mean, maintaining a certain rms-scatter which does not re- 
duce with time, i.e. an efficient relaxation to the medium 
temperature docs not occur. Note that the absolute size 
of the scatter is larger in the hot part of the slab. This 
is a result of the timestep criterion (25), which manages to 
hedge the rms-noise to something of order ~ aT. Reducing 
the timestep parameter a can thus improve the behaviour, 
while simultaneously increasing the computational cost sig- 
nificantly. For coarse timestepping (or high conductivities) 
the integration with this method can easily become unsta- 
ble. 



a deviation of an individual particle's temperature from the 
local mean is decaying only very slowly. 

The mixed formulation (C), shown in the bottom right, 
obviously shows the best behaviour in this test. It suppresses 
noise quickly, matches the analytical solution very well, and 
allows the largest timesteps of all schemes we tested. In 
this formulation, particles which have a large deviation from 
the local average temperature try to equalise this difference 
quickly, while a particle that is already close to the mean 
is not 'pulled away' by neighbours that have largo devia- 
tions. Apparently, this leads to better behaviour than for 
schemes (A) and (B), particularly when individual timesteps 
are used. We hence choose formulation (C) as our default 
method. 



Formulation (B), which uses the smoothed kernel- 
interpolated temperature field for both particles in each 
pair, docs significantly better in this respect (bottom left 
panel). However, the particle temperatures only very slowly 
approach the local mean value and hence the analytical solu- 
tion in this case. This is because the smoothing here is quite 
efficient in eliminating the small-scale noise, meaning that 



3.2 Conduction in a three-dimensional problem 

As a simple test of conduction in an intrinsically three- 
dimensional problem, we consider the temporal evolution 
of a point-like thermal energy perturbation. The time evolu- 
tion of an initial (5-function is given by the three-dimensional 
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Figure 3. Time evolution of the temperature field in an elementary three-dimensional conduction problem, displayed at two different 
times. We here consider the spreading of a narrow Gaussian temperature profile, which corresponds to the Green's function of the 

conduction problem for constant conductivity. Dots show the specific energies of individual simulation ijarticles, while the solid line 
marks the analytical result. Relative deviations arc shown in the lower panels. The numerical result maintains the Gaussian profile very 
well at all times. At early times, when the profile is sampled with few particles, a small reduction in the effective speed of conduction is 
seen, which however becomes increasingly unimportant at later times. 



Green's function for the conduction problem, 

where Cv = u/T is the heat capacity. We again consider 
a solid material, realized with a glass configuration of 50^ 
simulation particles with a mass of 1 g each and a moan par- 
ticle spacing of 1 cm. We choose a conductivity of kc„ = 
Icm^s"^, as before. We give the material a specific energy 
per unit mass of mo ~ 1000 erg g~^, and add a perturbation 
of 10000.0 erg cm^ g"^ x G{x, y, z, to = 10 sec). These initial 
conditions correspond to a (5-function perturbation that has 
already evolved for a brief period; this should largely elimi- 
nate timing offsets that would arise in the later evolution if 
we let the intrinsic SPH smoothing wash out a true initial 
(5-function. 

Figure 3 shows the specific energy profile after evolving 
the conduction problem for an additional 10 or 20 seconds, 
respectively. In the two bottom panels, we show the rela- 
tive deviation from the analytic solution. Given that quite 
coarse timestepping was used for this problem (about ~ 32 
timesteps for 20 sec), the match between the theoretical and 
numerical solutions is quite good. Even at a radial distance 
of 20 cm, where the Gaussian profile has dropped to less 
than one hundredth of its central amplitude, relative devi- 
ations are at around 10 percent for t — 20 sec and below 3 
percent at t = 30 sec. Note however that at all times the 
numerical solution maintains a nice Gaussian shape so that 
the deviations can be interpreted as a small modification in 
the effective conductivity. We then see that at late times, 
when the temperature gradients are resolved by more par- 
ticles, the SPH estimate of the conductivity becomes ever 
more accurate. 



4 SPHERICAL MODELS FOR CLUSTERS OF 
GALAXIES 

Observed temperature profiles of clusters of galaxies are of- 
ten characterised by a central decline of temperature, while 
otherwise appearing fairly isothermal over the measurable 
radial range. Provided the conductivity is non-negligible, 
there should hence be a conductive heat flux into the inner 
parts, which would then counteract central cooling losses. 
Motivated by this observation, Zakamska & Narayan (2003) 
(ZN henceforth) have constructed simple analytic models for 
the structure of clusters, invoking as a key assumption a lo- 
cal equilibrium between conduction and cooling. Combined 
with the assumption of hydrostatic equilibrium and spheri- 
cal symmetry, they were able to quite successfully reproduce 
the temperature profiles of a number of observed clusters. 

We here use the model of ZN as a test-bed to check the 
validity of our conduction modelling in a realistic cosmo- 
logical situation. In addition, the question for how long the 
ZN solution can be maintained is of immediate interest, and 
we address this with our simulations of spherical clusters as 
well. 

For definiteness, we briefly summarise the method of ZN 
for constructing a cluster equilibrium model. The cluster is 
assumed to be spherically symmetric, with the g£is of density 
p(r) and pressure P{r) being in hydrostatic equilibrium, i.e. 

1 dP d* 

-X- = -l^' (28) 
p ar ar 

where the gravitational potential $ is given by 

The dark matter density pdm(t') is described by a standard 
NFW (Navarro, Frenk & White, 1996) halo, optionally mod- 
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Figure 4. Local cooling (solid line) and conductive heating rates 
(dotted line) as a function of radius in our model for the cluster 
A2390, after a time of t = 0.1-5 Gyr. It is seen that both energy 
transfer rates cancel each other with reasonable accuracy. 



ified by ZN with the introduction of a small softened dark 
matter core. The structure of the dark halo can then be fully 
specified by the virial radius r2oo, defined as the radius en- 
closing a mean overdensity of 200 with respect to the critical 
density, and a scale length r^. 

The radial heat flux F due to electron conduction is 
given by 



F = -K-T-, 
ar 



(30) 



where the conductivity k is taJten to be a constant fraction 
of the Spitzer value. In order to balance local cooling losses, 
we require 

1 d 



dr 



{rF) = -j, 



(31) 



where j is describing the local gas cooling rate. If cooling 
is dominated by thermal bremsstrahlung, j can be approxi- 
mated as 



j = 2.1 X 10 



1/2 



ergs cm s 



(32) 



The electron number density herein correlates with the 
gas density p and depends on the hydrogen mass fraction X 
we assume for the primordial matter. With a proton mass 
of rup, it equals rie = p{X + l)/(2mp). Following ZN, we use 
this simplification in setting up our cluster models. 

Equations (28) to (32) form a system of differential 
equations that can be integrated from inside out once ap- 
propriate boundary conditions are specified. Wo adopt the 
values for central gas density and temperature, total dark 
matter mass, and scale radius determined by ZN for the 
cluster A2390 such that the resulting temperature profile 
matches the observed one well in the range < r < Tiout, 
where i?out was taken to be 2 . 

Integrating the system of equations out to JZout, we ob- 
tain a result that very well matches that reported by ZN. 
However, we also need to specify the structure of the cluster 
in its outer parts in order to be able to simulate it as an iso- 
lated system. ZN simply assumed this part to be isothermal 
at the temperature Tout = T'(i?out), thereby implicitly as- 



suming that the equilibrium condition between cooling and 
conduction invoked for the inner parts is not valid any more. 
It is a bit unclear why such a sudden transition should occur, 
but the alternative assumption, that equation (31) holds out 
to the virial radius, clearly leads to an unrealistic global tem- 
perature profile. In this case, the temperature would have to 
monotonically increase out to the outermost radius. Also, 
since the cumulative radiative losses out to a radius r have 
to be balanced by a conductive heat flux of equal magnitude 
at that radius, the heat flux also monotonically increases, 
such that all the energy radiated by the cluster would have 
to be supplied to the cluster at the virial radius. 

We nevertheless examined both approaches, i.e. we con- 
struct cluster models where we solve equations (28)- (32) for 
the whole cluster out to a radius of r2oo, and secondly, we 
construct initial conditions where we follow ZN by truncat- 
ing the solution at 2rs, continuing into the outer parts with 
an isothermal solution that is obtained by dropping equa- 
tions (30) and (31). Since conduction may still be important 
at radii somewhat larger than 2rs, these two models may 
hence be viewed as bracketing the expected real behaviour 
of the cluster in a model where cooling is balanced by con- 
duction. 

For the plasma conductivity, we assumed a value of 
0.3 Ksp, as proposed by ZN in their best-fit solution for 
A2390. Note that this value implies that a certain degree 
of suppression of conduction by magnetic fields is present, 
but that this effect is (perhaps optimistically) weak, corre- 
sponding to what is expected for chaotically tangled mag- 
netic fields (Narayan & Medvedev, 2001). 

Having obtained a solution for the static cluster model, 
we realized it as 3D initial conditions for GADGET. We used 
2 X 10^ gas particles with a total baryonic mass of around 
8.6 X 10^* Mq. For simplicity, we described the NFW dark 
matter halo of mass 2.4 x 10^^ Mq as a static potential. 
We then simulated the evolution of the gas subject to self- 
gravity (with a gravitational softening length of 3kpc), ra- 
diative cooling, and thermal conduction, but without allow- 
ing for star formation and associated feedback processes. 

By construction, we expect that thermal conduction will 
bo able to offset radiative cooling for those initial conditions, 
at least for some time. This can be verified in Figure 4, 
where we plot the local cooling rate and conductive heat- 
ing rate as a function of radius. Indeed, at time 0.15 Gyr, 
after a brief initial relaxation period, the two energy trans- 
fer rates exhibit the same magnitude and cancel each other 
with good accuracy. This represents a nice validation of our 
numerical implementation of thermal conduction with the 
temperature-dependent Spitzer-rate. 

As the simulation continues, we see that the core tem- 
perature of the cluster slowly drifts to somewhat higher tem- 
perature. A secular evolution of some kind is probably un- 
avoidable, since the balance between cooling and conductive 
heating cannot be perfectly static. The cluster of course still 
loses all the energy it radiates, which eventually must give 
rise to a slow quasi-static inflow of gas, and a corresponding 
change of the inner structure of the cluster. Because of the 
different temperature dependences of cooling and conductive 
heating, it is also not clear that the balance between cool- 
ing and conductive heating represents a stable dynamical 
state. While a stability analysis by ZN and Kim & Narayan 
(2003) suggests that thermal instability is sufficiently sup- 
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Figure 5. Temperature profile (left) and cumulative mass profile (right) of our model for the cluster A2390, after a simulated time 

of 0.6 Gyr. We compare simulations with (triangles) and without (diamonds) thermal conduction using a conductivity of O.Sftsp- It is 
seen tiiat without conduction, the inner regions of the cluster cool down significantly, causing mass to sink towards the centre as central 
pressure support is partially lost. If thermal conduction is included, the central cooling losses arc offset by conductive heating from outer 
regions, preventing any significant change of the baryonic mass profile. In fact, in this case we even observe a slow secular evolution 
towards higher temperature in the inner parts. 



pressed in models with conduction, Soker (2003) argues that 
local perturbations in the cooling flow region will grow to 
the non-linear regime rather quickly and that a steady so- 
lution with a constant heat conduction may therefore not 
exist. 

In any case, it is clear that the inclusion of conduction 
strongly reduces central mass drop out due to cooling. This 
is demonstrated explicitly in Figure 5, where we show the 
temperature and cumulative baryonic mass profiles of the 
cluster after a time of 0.6 Gyr. We compare these profiles 
also to an identical simulation where conduction was not 
included. Unlike the conduction run, this cooling-only sim- 
ulation shows a very substantial modification of its temper- 
ature profile in the inner parts, where the core region inside 
~ 0.1 Mpc becomes much colder than the initial state. We 
can also see that the cooling-only simulation has seen sub- 
stantial baryonic inflow as a result of central mass drop-out. 
Within 0.6 Gyr, the baryonic mass enclosed in a sphere with 
radius 50 kpc has increased by 86 percent in this simulation, 
clearly forming a cooling flow. This contrasts strongly with 
the run that includes conduction, where beyond a cental dis- 
tance of ~ 0.02 Mpc, there is no significant difference in the 
cumulative mass profile between t = and t = 0.6 Gyr. 

We also checked that runs with a reduced conductivity 
show a likewise reduced effect, suppressing the core cooling 
and matter inflow to a lesser extent than with our default 
model with 0.3 Ksp. 

Finally, we note that the results discussed here are quite 
insensitive to whether we set-up the initial conditions with 
an isothermal outer part, or whether we continue the equi- 
librium solution to the virial radius. There is only a small 
difl^erence in the long term evolution of the solution. 

The model presented here helps us to verify the robust- 
ness of our conduction implementation when a temperature 
dependent conductivity is used, but it does not account for 
the time dependence of conduction expected in the cosmo- 
logical context. The static, spherically symmetric construc- 
tion does not reflect the hierarchical growth of clusters that 



is a central element of currently favoured cosmologies. To be 
able to make reliable statements about the effects of ther- 
mal conduction on galaxy clusters, one therefore has to trace 
their evolution from high redshift to the present, which is 
best done in full cosmological simulations. Here, the hier- 
archical growth of clusters implies that the conductivity is 
becoming low early on, when the ICM temperature is low, 
and becomes only large at late times when the cluster forms. 



5 COSMOLOGICAL CLUSTER SIMULATIONS 

In this section, we apply our new numerical scheme for the 
treatment of thermal conduction in fully self-consistent cos- 
mological simulations of cluster formation. We will in par- 
ticular compare the results of simulations with and without 
conduction, both for runs that follow only adiabatic gas- 
dynamics, and runs that also include radiative cooling of 
gas. While it is beyond the scope of this work to present a 
comprehensive analysis of the efi^ects of conduction in cos- 
mological simulations, wc here want to investigate a sot of 
small fiducial runs in order to further validate the stability 
of our method for real-world applications, and secondly, to 
give a first fiavour of the expected effects of conduction in 
simulated clusters. A more detailed analysis of cosmologi- 
cal implications of conduction is presented in a companion 
paper (Dolag et al., 2004). 

We focus on a single cluster, extracted from the 'GIF' 
simulation (Kauffmann et al., 1999) of the ACDM model, 
and resimulated using the "zoomed initial conditions" tech- 
nique (Tormen et al., 1997). To this end, the particles that 
make up the cluster at the present time are traced back 
to their original coordinates in the initial conditions. The 
Lagrangiaii region of the cluster identified in this way is 
then resampled with particles of smaller mass, and addi- 
tional small-scale perturbations from the CDM power spec- 
trum are added appropriately. Far away from the cluster. 
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Figure 6. Projections of mass-weighted temperature (left column), X-ray emissivity (middle column), and gas mass density (right 
column) for our cluster simulations at 2 = 0.13. From the top to the bottom row, we show the same cluster but simulated with different 
physics: Adiabatic gasdynamics only, adiabatic plus thermal conduction, radiative cooling and star formation, and finally, cooling, star 
formation and conduction. Each panel displays the gas contained in a box of side-length 8.6 Mpc centred on the cluster. Full Spitzer 
conductivity was assumed. 
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Figure 7. Cumulative baryonic mass profile (top), temperature 
profile (middle) and X-ray cmissivity profile (bottom) of the sim- 
ulated cluster at z = 0.13. In each panel, we compare the same 
cluster simulation run with different physical models for the gas: 
Adiabatic gasdynamics only, adiabatic plus thermal conduction, 
radiative cooling and star formation without conduction, and fi- 
nally, cooling, star formation and conduction. For the models in- 
cluding conduction, full Spitzer conductivity was assumed. Note 
that the X-ray emissivity is plotted such that the area under the 
curve is proportional to the total bolometric X-ray luminosity. 



the resolution is progressively degraded by using particles of 
ever larger mass. 

The cluster we selected has a virial mass of 1.1 x 
10^^ M0. It is the same cluster considered in the high- 
resolution study of Springel et al. (2001a), with our mass 
resolution corresponding to their SI simulation, except that 
we split each of the ~ 450000 high resolution dark matter 
particles into a gas and a dark matter particle (assuming 
fib = 0.04), yielding a gas mass resolution of 8.4 x 10^ M©. 
The boundary region was sampled with an additional 3 mil- 
lion dark matter particles. We then evolved the cluster for- 
ward in time from a starting redshift z = 30 to the present, 
z = Q, using a comoving gravitational softening length of 
20 kpc. 

We consider 4 different simulations of the cluster, each 
with different physical models for the gas: (1) Adiabatic gas- 
dynamics only, (2) adiabatic gas and conduction, (3) radia- 
tive cooling and star formation without conduction, and fi- 
nally, (4) radiative cooling, star formation and conduction. 
The simulations with cooling and star formation use the 
sub-resolution model of Springel & Hcrnquist (2003) for the 
multi-phase structure of the ISM (without including the op- 
tional feedback by galactic winds offered by the model). In 
the two runs with conduction, wo adopt the full Spitzer rate 
for the conductivity. While this value is unrealistically large 
for magnetised clusters, it serves our purpose here in high- 
lighting the effects of conduction, given also that our cluster 
is not particularly hot, so that effects of conductivity can be 
expected to be weaker than in very rich clusters. 

In Figure 6, we show projections of the mass-weighted 
temperature. X-ray emissivity, and gas mass density for all 
four simulations. Each panel displays the gas contained in a 
box of side-length 8.6 Mpc centred on the cluster. Compar- 
ing the simulations with and without conduction, it is nicely 
seen how conduction tends to wipe out small-scale temper- 
ature fluctuations. It is also seen that the outer parts of the 
cluster become hotter when conduction is included. 

These trends axe also borne out quantitatively when 
studying radial profiles of cluster properties in more detail. 
In Figure 7, wo compare the cumulative baryonic mass pro- 
file, the temperature profile, and the radial profile of X-ray 
emission for all four simulations. Note that the innermost 
bins, for R < 30 kpc, may be affected by numerical reso- 
lution effects. Interestingly, the temperature profiles of the 
runs with conduction are close to being perfectly isother- 
mal in the inner parts of the cluster. While this does not 
represent a large change for the adiabatic simulation, which 
is close to isothermal anyway, the simulation with radia- 
tive cooling is changed significantly. Without conduction, 
the radiative run actually shows a pronounced rise in the 
temperature profile in the range of 100 — 200 kpc, as a result 
of compressional heating when gets flows in to replace gas 
that is cooling out of the ICM in a cooling flow. Only in the 
innermost regions, where cooling becomes rapid, we see a 
distinct drop of the temperature. Interestingly, conduction 
eliminates this feature in the temperature profile, by trans- 
porting the corresponding heat energy from the max;imum 
both to parts of the cluster further out and to the innermost 
parts. The latter effect is probably small, however, because a 
smooth decline in the temperature profile in the inner parts 
of the cluster, as present in the ZN model, does not appear 
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in the simulation. As a consequence, a strong conductive 
heat flow from outside to inside cannot develop. 

Conduction may also induce changes in the X-ray emis- 
sion of the clusters, which we show in the bottom panel of 
Figure 7. Interestingly, the inclusion of conduction in the 
adiabatic simulation has a negligible effect on the X-ray lu- 
minosity. This is because in contrast to previous suggestions 
(Loeb, 2002), the cluster does not lose a significant fraction 
of its thermal energy content to the outside intergalactic 
medium, and the changes in the relevant part of the gas and 
temperature profile are rather modest. We do note however 
that the redistribution of thermal energy within the clus- 
ter leads to a substantial increase of the temperature of the 
outer parts of the cluster. 

For the simulations with cooling, the changes of the 
X-ray properties are more significant. Interestingly, we find 
that allowing for thermal conduction leads to a net increase 
of the bolo metric luminosity of our simulated cluster. The 
panel with the cumulative baryon mass profile reveals that 
conduction is also ineffective in significantly suppressing the 
condensation of mass in the core regions of the cluster. In 
fact, it may even lead to the opposite effect. We think this 
behaviour simply occurs because a temperature profile with 
a smooth decline towards the centre, which would allow the 
conductive heating of this part of the cluster, is not forming 
in the simulation. Instead, the conductive heat flux is point- 
ing primarily from the inside to the outside, which may then 
be viewed as an additional "cooling" process for the inner 
cluster regions. 

It is interesting to note that in spite of the structural ef- 
fects that thermal conduction has on the ICM of the cluster, 
it does not affect its star formation history significantly. In 
the two simulations that include cooling and star formation, 
the stellar component of the mass profile (Figure 7) does 
not show any sizeable difference between the run including 
thermal conduction and the one without it. This does not 
come as surprise as 90% of the stellar content of the clus- 
ter has formed before a rcdshift of z — 0.85. At these early 
times, the temperature of the gas in the protocluster was 
much lower, such that conduction was unimportant. In fact, 
for that reason, the stellar mass profiles for both simulation 
runs coincide. 

In summary, our initial results for this cluster suggest 
that conduction can be important for the ICM, provided 
the effective conductivity is a sizable fraction of the Spitzer 
value. However, the interplay between radiative cooling and 
conduction is clearly complex, and it is presently unclear 
whether temperature profiles like those observed can arise 
in self-consistent cosmological simulations. We caution that 
one should not infer too much from the single object we 
examined here. A much larger set of cluster simulations will 
be required to understand this topic better. 



moderating cooling fiows in clusters of galaxies. Such sim- 
ulations are then an ideal tool to make reliable predictions 
of the complex interplay between the nonlinear processes of 
cooling and conduction during structure formation. 

In this paper, wc have presented a detailed numerical 
methodology for the treatment of conduction in cosmological 
SPH simulations. By construction, our method manifestly 
conserves thermal energy, and we have formulated it such 
that it is robust against the presence of small-scale temper- 
ature noise. Wc have implemented this method in a modern 
parallel code, capable of carrying out large, high-resolution 
cosmological simulations. 

Using various test problems, wc have demonstrated the 
accuracy and robustness of our rmmerical scheme for con- 
duction. We then applied our code to a first set of cosmo- 
logical cluster formation simulations, comparing in particu- 
lar simulations with and without conduction. While these 
results are preliminary, they already hint that the phe- 
nomenology of the coupled dynamics of radiative cooling 
and conduction is complex, and may give rise to results that 
were perhaps not anticipated by earlier analytic modelling 
of static cluster configurations. 

For example, wc found that conduction docs not neces- 
sarily reduce a central cooling flow in our simulations; the 
required smoothly declining temperature profile in the in- 
ner cluster regions does not readily form in our cosmological 
simulations. Instead, the profiles we find are either flat, or 
have a tendency to slightly rise towards the centre, akin to 
what is seen in the cooling-only simulations. In this situa- 
tion, conduction may in fact lead to additional cooling in 
certain situations, by either transporting thermal energy to 
the outer parts, or by modifying the temperature and den- 
sity structure in the relevant parts of the cluster such that 
cooling is enhanced. This can then manifest itself in an in- 
crease of the bolometric X-ray luminosity at certain times, 
which is actually the case for our model cluster at z = 0. 
Interestingly, we do not find that our cluster loses a signifi- 
cant fraction of its thermal heat content by conducting it to 
the external intergalactic medium. 

In a companion paper (Dolag et al., 2004), we analyse a 
larger set of cosmological cluster simulations, computed with 
much higher resolution and with more realistic sub-Spitzer 
conductivities. This set of cluster allows us to investigate, 
e.g., conduction effects as a function of cluster temperature 
and the influence of conduction on cluster scaling relations. 
While our first results of this work suggest that conduc- 
tion by itself may not resolve the cooling-flow puzzle, it also 
shows that conduction has a very strong influence on the 
thermodynamic state of rich clusters if the effective conduc- 
tivity is a small fraction of the Spitzer value or more. In 
future work, it will hence be very interesting and impor- 
tant to understand the rich phenomenology of conduction 
in clusters in more detail. 



6 CONCLUSIONS 

Hot plasmas like those found in clusters of galaxies are effi- 
ciently conducting heat, unless electron thermal conduction 
is heavily suppressed by magnetic fields. Provided the lat- 
ter is not the case, heat conduction should therefore be in- 
cluded in hydrodynamical cosmological simulations, given, 
in particular, that conduction could play a decisive role in 
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